Data analysis for:

Why are telomeres the length that they are? Insight from a phylogenetic comparative analysis

Workflow of the analyses produced by Dylan J. Padilla Perez, School of Life Sciences, Arizona State University, Tempe, AZ 85287, USA.

Library

library(AICcmodavg)
library(ape)
library(caper)
library(car)
library(coda)
library(extrafont)
library(geiger)
library(kableExtra)
library(MuMIn)
library(nlme)
library(pbapply)
library(phylopath)
library(phytools)
library(plotrix)
library(rphylopic)
library(scales)
library(shape)
library(xtable)
R.version
>                _                           
> platform       x86_64-apple-darwin17.0     
> arch           x86_64                      
> os             darwin17.0                  
> system         x86_64, darwin17.0          
> status                                     
> major          4                           
> minor          2.2                         
> year           2022                        
> month          10                          
> day            31                          
> svn rev        83211                       
> language       R                           
> version.string R version 4.2.2 (2022-10-31)
> nickname       Innocent and Trusting
sessionInfo()
> R version 4.2.2 (2022-10-31)
> Platform: x86_64-apple-darwin17.0 (64-bit)
> Running under: macOS Catalina 10.15.7
> 
> Matrix products: default
> BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
> LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
> 
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
> 
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base     
> 
> other attached packages:
>  [1] rmarkdown_2.21   xtable_1.8-4     shape_1.4.6      scales_1.2.1    
>  [5] rphylopic_1.2.0  plotrix_3.8-2    phylopath_1.1.3  pbapply_1.7-2   
>  [9] nlme_3.1-162     MuMIn_1.47.5     kableExtra_1.3.4 geiger_2.0.11   
> [13] phytools_2.0-3   maps_3.4.1       extrafont_0.19   coda_0.19-4     
> [17] car_3.1-2        carData_3.0-5    caper_1.0.1      mvtnorm_1.1-3   
> [21] MASS_7.3-59      ape_5.7-1        AICcmodavg_2.3-2
> 
> loaded via a namespace (and not attached):
>  [1] subplex_1.8             doParallel_1.0.17       webshot_0.5.4          
>  [4] httr_1.4.5              numDeriv_2016.8-1.1     tools_4.2.2            
>  [7] bslib_0.4.2             utf8_1.2.3              R6_2.5.1               
> [10] colorspace_2.1-0        tidyselect_1.2.0        mnormt_2.1.1           
> [13] phangorn_2.11.1         grImport2_0.2-0         curl_5.0.0             
> [16] compiler_4.2.2          extrafontdb_1.0         cli_3.6.1              
> [19] rvest_1.0.3             expm_0.999-7            xml2_1.3.4             
> [22] unmarked_1.2.5          sass_0.4.5              quadprog_1.5-8         
> [25] systemfonts_1.0.4       stringr_1.5.0           digest_0.6.31          
> [28] svglite_2.1.1           base64enc_0.1-3         jpeg_0.1-10            
> [31] pkgconfig_2.0.3         htmltools_0.5.5         fastmap_1.1.1          
> [34] highr_0.10              rlang_1.1.1             rstudioapi_0.14        
> [37] VGAM_1.1-8              optimParallel_1.0-2     farver_2.1.1           
> [40] jquerylib_0.1.4         generics_0.1.3          combinat_0.0-8         
> [43] jsonlite_1.8.4          dplyr_1.1.2             magrittr_2.0.3         
> [46] Matrix_1.5-4            Rcpp_1.0.10             munsell_0.5.0          
> [49] fansi_1.0.4             abind_1.4-5             lifecycle_1.0.3        
> [52] scatterplot3d_0.3-43    stringi_1.7.12          yaml_2.3.7             
> [55] clusterGeneration_1.3.7 plyr_1.8.8              grid_4.2.2             
> [58] parallel_4.2.2          lattice_0.21-8          splines_4.2.2          
> [61] knitr_1.42              pillar_1.9.0            igraph_1.4.2           
> [64] codetools_0.2-19        stats4_4.2.2            fastmatch_1.1-3        
> [67] XML_3.99-0.14           glue_1.6.2              evaluate_0.20          
> [70] deSolve_1.35            vctrs_0.6.2             png_0.1-8              
> [73] foreach_1.5.2           Rttf2pt1_1.3.12         gtable_0.3.3           
> [76] cachem_1.0.8            ggplot2_3.4.2           xfun_0.39              
> [79] rsvg_2.4.0              survival_3.5-5          viridisLite_0.4.2      
> [82] tibble_3.2.1            iterators_1.0.14
set.seed(80)

## Dataset

data <- read.csv("new_data_2024.csv")
str(data)
> 'data.frame': 159 obs. of  18 variables:
>  $ Common.name               : chr  "Adelie penguin" "Alpine swift " "American bald eagle " "American flamengo" ...
>  $ Scientific_name           : chr  "Pygoscelis_adeliae" "Tachymarptis_melba" "Haliaeetus_leucocephalus" "Phoenicopterus_ruber" ...
>  $ Domain                    : chr  "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
>  $ Kingdom                   : chr  "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
>  $ Phylum                    : chr  "Chordata" "Chordata" "Chordata" "Chordata" ...
>  $ Class                     : chr  "Aves" "Aves" "Aves" "Aves" ...
>  $ Order                     : chr  "Sphenisciformes" "Apodiformes" "Accipitriformes" "Phoenicopteriformes" ...
>  $ Family                    : chr  "Spheniscidae" "Apodidae" "Accipitridae" "Phoenicopteridae" ...
>  $ Genus                     : chr  "Pygoscelis" "Tachymarptis" "Haliaeetus" "Phoenicopterus" ...
>  $ Species                   : chr  "adeliae" "melba" "leucocephalus" "ruber" ...
>  $ Endo_ectotherm            : chr  "endo" "endo" "endo" "endo" ...
>  $ Adult_mass_grams          : num  5020 103 3175 3066 117 ...
>  $ log_mass                  : num  3.7 2.01 3.5 3.49 2.07 ...
>  $ Lifespan_years            : num  20 26 48 44 17 25 16 28.2 29.9 43.7 ...
>  $ Loglifespan               : num  1.3 1.41 1.68 1.64 1.23 ...
>  $ Average_Telomere_Length_kb: num  7.35 18.78 2.1 18.31 30.2 ...
>  $ logTL                     : num  0.866 1.274 0.322 1.263 1.48 ...
>  $ Telomerase_activity       : int  NA NA NA NA NA NA NA NA NA NA ...
head(data)
>            Common.name          Scientific_name    Domain Kingdom   Phylum
> 1       Adelie penguin       Pygoscelis_adeliae Eukaryota Metazoa Chordata
> 2        Alpine swift        Tachymarptis_melba Eukaryota Metazoa Chordata
> 3 American bald eagle  Haliaeetus_leucocephalus Eukaryota Metazoa Chordata
> 4    American flamengo     Phoenicopterus_ruber Eukaryota Metazoa Chordata
> 5   American kestrel           Falco_sparverius Eukaryota Metazoa Chordata
> 6       Andouin's gull    Ichthyaetus_audouinii Eukaryota Metazoa Chordata
>   Class               Order           Family          Genus       Species
> 1  Aves     Sphenisciformes     Spheniscidae     Pygoscelis       adeliae
> 2  Aves         Apodiformes         Apodidae   Tachymarptis         melba
> 3  Aves     Accipitriformes     Accipitridae     Haliaeetus leucocephalus
> 4  Aves Phoenicopteriformes Phoenicopteridae Phoenicopterus         ruber
> 5  Aves       Falconiformes       Falconidae          Falco    sparverius
> 6  Aves     Charadriiformes          Laridae    Ichthyaetus     audouinii
>   Endo_ectotherm Adult_mass_grams log_mass Lifespan_years Loglifespan
> 1           endo           5020.0 3.700704             20    1.301030
> 2           endo            102.7 2.011570             26    1.414973
> 3           endo           3175.0 3.501744             48    1.681241
> 4           endo           3066.0 3.486572             44    1.643453
> 5           endo            117.0 2.068186             17    1.230449
> 6           endo            535.0 2.728354             25    1.397940
>   Average_Telomere_Length_kb     logTL Telomerase_activity
> 1                      7.350 0.8662873                  NA
> 2                     18.783 1.2737650                  NA
> 3                      2.100 0.3222193                  NA
> 4                     18.310 1.2626883                  NA
> 5                     30.200 1.4800069                  NA
> 6                     26.550 1.4240645                  NA
unique(data$Class)
> [1] "Aves"     "Fish"     "Mammalia" "Reptilia"
dat <- data
names(dat)
>  [1] "Common.name"                "Scientific_name"           
>  [3] "Domain"                     "Kingdom"                   
>  [5] "Phylum"                     "Class"                     
>  [7] "Order"                      "Family"                    
>  [9] "Genus"                      "Species"                   
> [11] "Endo_ectotherm"             "Adult_mass_grams"          
> [13] "log_mass"                   "Lifespan_years"            
> [15] "Loglifespan"                "Average_Telomere_Length_kb"
> [17] "logTL"                      "Telomerase_activity"
dat$log_mass <- log1p(dat$Adult_mass_grams)
dat$log_mass
>   [1]  8.5213844  4.6415021  8.0633778  8.0284552  4.7706846  6.2841342
>   [7]  2.9601051  7.4312997  6.0602907  8.0811658  5.7620514  4.6463121
>  [13]  6.6603190  4.7957905  3.9796817  6.1758673  7.4809922  4.3782696
>  [19]  6.4393504  2.5649494  7.2449415  2.9391619  9.0669318  8.9142228
>  [25]  5.5093883  7.4679423  3.8199077  6.6427475  2.7725887  5.2729996
>  [31]  6.8123451  8.3696208  6.9517722  5.9162021  7.2174434  2.8903718
>  [37]  6.9255952  4.9199809  2.6173958  3.0540012  8.2158178  6.8892246
>  [43]  6.8721281  4.8751973  2.9957323  8.8604992  2.5877640  8.7079788
>  [49]  2.5649494  6.2841342  6.3985949  2.8273136  8.3022658  0.5306283
>  [55]  8.0067008  7.7836406  7.3783837  9.2687036  0.4054651  1.0296194
>  [61]  1.5475625  1.0986123  0.8329091  2.5877640 10.8742854  8.7949764
>  [67]  8.0712185  6.5191473  9.4880479  8.4734503  9.4121377  8.0394799
>  [73]  8.0067008  8.2942996  7.2717037  8.5183925  7.4091364 10.5966597
>  [79]  5.4205350 18.4206808 17.1654147 13.5923683 11.2897944 11.5253579
>  [85] 11.0186455 12.2060776 12.9808021 11.7752974  9.7981826  9.9523253
>  [91] 13.5278298 12.2783980 12.1007177 11.0509059 10.5966597 13.0710722
>  [97]  8.2689882  8.3723986 11.6927522  9.4727816  2.3978953  2.6026897
> [103]  5.8607862  6.7345917  8.0067008  6.2166061  1.9878743  2.3025851
> [109]  8.1889669  8.3371091  7.9013774  7.4960973  4.6151205  3.7135721
> [115]  4.0775374 14.5925397 12.8584004 12.6115411 12.4292202 10.2576945
> [121]  8.7404967 11.0354701  8.8913740  6.8308742 10.5947830 11.0740483
> [127] 10.7140844 11.8482756  7.8461988  8.9763886  9.0162701  8.7584124
> [133] 15.3841267 14.9717629  9.9159595 10.9151066  6.2803958  7.0264268
> [139]  3.0680529  3.5835189  5.7071103  3.0680529  5.2922993  5.3033049
> [145]  9.5888453 11.9183972  1.7917595  3.0445224  9.2679487  5.7868974
> [151]  7.2647302  5.9210421  6.4669220  8.9683962  3.9926809  4.6634391
> [157]  7.2174434  8.2942996  4.5747110
dat[dat$Average_Telomere_Length_kb > 200, ]
>       Common.name Scientific_name    Domain Kingdom   Phylum    Class
> 107 Iberian shrew Sorex_granarius Eukaryota Metazoa Chordata Mammalia
>            Order    Family Genus   Species Endo_ectotherm Adult_mass_grams
> 107 Eulipotyphla Soricidae Sorex granarius           endo              6.3
>     log_mass Lifespan_years Loglifespan Average_Telomere_Length_kb   logTL
> 107 1.987874              3   0.4771213                        213 2.32838
>     Telomerase_activity
> 107                   1
#str(dat)

dat <- dat[!dat$Scientific_name == "Sorex_granarius", ] ## potential outlier


rownames(dat) <- dat$Scientific_name
head(dat)
>                                   Common.name          Scientific_name
> Pygoscelis_adeliae             Adelie penguin       Pygoscelis_adeliae
> Tachymarptis_melba              Alpine swift        Tachymarptis_melba
> Haliaeetus_leucocephalus American bald eagle  Haliaeetus_leucocephalus
> Phoenicopterus_ruber        American flamengo     Phoenicopterus_ruber
> Falco_sparverius           American kestrel           Falco_sparverius
> Ichthyaetus_audouinii          Andouin's gull    Ichthyaetus_audouinii
>                             Domain Kingdom   Phylum Class               Order
> Pygoscelis_adeliae       Eukaryota Metazoa Chordata  Aves     Sphenisciformes
> Tachymarptis_melba       Eukaryota Metazoa Chordata  Aves         Apodiformes
> Haliaeetus_leucocephalus Eukaryota Metazoa Chordata  Aves     Accipitriformes
> Phoenicopterus_ruber     Eukaryota Metazoa Chordata  Aves Phoenicopteriformes
> Falco_sparverius         Eukaryota Metazoa Chordata  Aves       Falconiformes
> Ichthyaetus_audouinii    Eukaryota Metazoa Chordata  Aves     Charadriiformes
>                                    Family          Genus       Species
> Pygoscelis_adeliae           Spheniscidae     Pygoscelis       adeliae
> Tachymarptis_melba               Apodidae   Tachymarptis         melba
> Haliaeetus_leucocephalus     Accipitridae     Haliaeetus leucocephalus
> Phoenicopterus_ruber     Phoenicopteridae Phoenicopterus         ruber
> Falco_sparverius               Falconidae          Falco    sparverius
> Ichthyaetus_audouinii             Laridae    Ichthyaetus     audouinii
>                          Endo_ectotherm Adult_mass_grams log_mass
> Pygoscelis_adeliae                 endo           5020.0 8.521384
> Tachymarptis_melba                 endo            102.7 4.641502
> Haliaeetus_leucocephalus           endo           3175.0 8.063378
> Phoenicopterus_ruber               endo           3066.0 8.028455
> Falco_sparverius                   endo            117.0 4.770685
> Ichthyaetus_audouinii              endo            535.0 6.284134
>                          Lifespan_years Loglifespan Average_Telomere_Length_kb
> Pygoscelis_adeliae                   20    1.301030                      7.350
> Tachymarptis_melba                   26    1.414973                     18.783
> Haliaeetus_leucocephalus             48    1.681241                      2.100
> Phoenicopterus_ruber                 44    1.643453                     18.310
> Falco_sparverius                     17    1.230449                     30.200
> Ichthyaetus_audouinii                25    1.397940                     26.550
>                              logTL Telomerase_activity
> Pygoscelis_adeliae       0.8662873                  NA
> Tachymarptis_melba       1.2737650                  NA
> Haliaeetus_leucocephalus 0.3222193                  NA
> Phoenicopterus_ruber     1.2626883                  NA
> Falco_sparverius         1.4800069                  NA
> Ichthyaetus_audouinii    1.4240645                  NA
## Trees

full_data_tree <- read.tree("spp.tree.nwk")
is.ultrametric(full_data_tree)
> [1] FALSE
full_data_tree <- force.ultrametric(full_data_tree)
> ***************************************************************
> *                          Note:                              *
> *    force.ultrametric does not include a formal method to    *
> *    ultrametricize a tree & should only be used to coerce    *
> *   a phylogeny that fails is.ultrametric due to rounding --  *
> *    not as a substitute for formal rate-smoothing methods.   *
> ***************************************************************
is.ultrametric(full_data_tree)
> [1] TRUE
full_data_tree
> 
> Phylogenetic tree with 158 tips and 157 internal nodes.
> 
> Tip labels:
>   Raja_montagui, Torpedo_marmorata, Urolophus_paucimaculatus, Squalus_acanthias, Mustelus_asterias, Callorhinchus_milii, ...
> Node labels:
>   , '36', '37', '13', '14', '25', ...
> 
> Rooted; includes branch lengths.
## Full_tree

full_data_tree$tip.label
>   [1] "Raja_montagui"              "Torpedo_marmorata"         
>   [3] "Urolophus_paucimaculatus"   "Squalus_acanthias"         
>   [5] "Mustelus_asterias"          "Callorhinchus_milii"       
>   [7] "Oncorhynchus_mykiss"        "Lutjanus_argentimaculatus" 
>   [9] "Acanthopagrus_schlegelii"   "Dicentrarchus_labrax"      
>  [11] "Macquaria_ambigua"          "Platycephalus_bassensis"   
>  [13] "Trachurus_novaezelandiae"   "Pseudocaranx_wrighti"      
>  [15] "Oreochromis_niloticus"      "Oryzias_latipes"           
>  [17] "Fundulus_heteroclitus"      "Xiphophorus_hellerii"      
>  [19] "Xiphophorus_couchianus"     "Xiphophorus_maculatus"     
>  [21] "Nothobranchius_kadleci"     "Nothobranchius_furzeri"    
>  [23] "Nothobranchius_rachovii"    "Upeneichthys_vlamingii"    
>  [25] "Gadus_morhua"               "Carassius_auratus"         
>  [27] "Cyprinus_carpio"            "Danio_rerio"               
>  [29] "Anguilla_rostrata"          "Didelphis_virginiana"      
>  [31] "Sorex_granarius"            "Sorex_araneus"             
>  [33] "Atelerix_albiventris"       "Crocuta_crocuta"           
>  [35] "Felis_catus"                "Panthera_tigris"           
>  [37] "Canis_lupus"                "Meles_meles"               
>  [39] "Ailurus_fulgens"            "Zalophus_californianus"    
>  [41] "Ursus_maritimus"            "Equus_caballus"            
>  [43] "Equus_grevyi"               "Tapirus_indicus"           
>  [45] "Ceratotherium_simum"        "Camelus_dromedarius"       
>  [47] "Tursiops_truncatus"         "Balaena_mysticetus"        
>  [49] "Eschrichtius_robustus"      "Hexaprotodon_liberiensis"  
>  [51] "Bos_taurus"                 "Capra_hircus"              
>  [53] "Ovis_aries"                 "Muntiacus_muntjak"         
>  [55] "Muntiacus_reevesi"          "Rangifer_tarandus"         
>  [57] "Giraffa_camelopardalis"     "Sus_scrofa"                
>  [59] "Pteropus_rodricensis"       "Tadarida_brasiliensis"     
>  [61] "Myotis_lucifugus"           "Sylvilagus_aquaticus"      
>  [63] "Oryctolagus_cuniculus"      "Lepus_californicus"        
>  [65] "Ochotona_princeps"          "Aplodontia_rufa"           
>  [67] "Sciurus_carolinensis"       "Peromyscus_eremicus"       
>  [69] "Cricetulus_griseus"         "Rattus_norvegicus"         
>  [71] "Mus_musculus"               "Mus_spretus"               
>  [73] "Castor_canadensis"          "Heterocephalus_glaber"     
>  [75] "Hydrochoerus_hydrochaeris"  "Lemur_catta"               
>  [77] "Saimiri_sciureus"           "Ateles_geoffroyi"          
>  [79] "Gorilla_gorilla"            "Homo_sapiens"              
>  [81] "Pan_paniscus"               "Pan_troglodytes"           
>  [83] "Pongo_pygmaeus"             "Macaca_mulatta"            
>  [85] "Macaca_fascicularis"        "Macaca_nemestrina"         
>  [87] "Tupaia_tana"                "Tupaia_belangeri"          
>  [89] "Procavia_capensis"          "Loxodonta_africana"        
>  [91] "Elephas_maximus"            "Setifer_setosus"           
>  [93] "Galegeeska_rufescens"       "Macroscelides_proboscideus"
>  [95] "Chaetophractus_vellerosus"  "Myrmecophaga_tridactyla"   
>  [97] "Choloepus_hoffmanni"        "Liasis_fuscus"             
>  [99] "Thamnophis_sirtalis"        "Lacerta_agilis"            
> [101] "Zootoca_vivipara"           "Alligator_mississippiensis"
> [103] "Alligator_sinensis"         "Tachymarptis_melba"        
> [105] "Sula_sula"                  "Fregata_minor"             
> [107] "Pygoscelis_adeliae"         "Macronectes_giganteus"     
> [109] "Macronectes_halli"          "Fulmarus_glacialis"        
> [111] "Calonectris_diomedea"       "Oceanodroma_leucorhoa"     
> [113] "Diomedea_exulans"           "Thalassarche_melanophris"  
> [115] "Antigone_antigone"          "Sterna_hirundo"            
> [117] "Ichthyaetus_audouinii"      "Larus_fuscus"              
> [119] "Rissa_tridactyla"           "Uria_lomvia"               
> [121] "Cepphus_grylle"             "Calidris_alpina"           
> [123] "Calidris_pugnax"            "Haematopus_ostralegus"     
> [125] "Phoenicopterus_ruber"       "Buteo_swainsoni"           
> [127] "Buteo_jamaicensis"          "Haliaeetus_leucocephalus"  
> [129] "Accipiter_gentilis"         "Gyps_fulvus"               
> [131] "Amazona_amazonica"          "Strigops_habroptila"       
> [133] "Acrocephalus_sechellensis"  "Hirundo_rustica"           
> [135] "Riparia_riparia"            "Tachycineta_albilinea"     
> [137] "Tachycineta_bicolor"        "Parus_major"               
> [139] "Turdus_merula"              "Taeniopygia_guttata"       
> [141] "Lonchura_striata"           "Chloebia_gouldiae"         
> [143] "Coloeus_monedula"           "Aphelocoma_ultramarina"    
> [145] "Aphelocoma_coerulescens"    "Passerculus_sandwichensis" 
> [147] "Vireo_olivaceus"            "Falco_sparverius"          
> [149] "Branta_leucopsis"           "Chrysolophus_pictus"       
> [151] "Meleagris_gallopavo"        "Gallus_gallus"             
> [153] "Colinus_virginianus"        "Trachemys_scripta"         
> [155] "Chrysemys_picta"            "Emys_orbicularis"          
> [157] "Chelonia_mydas"             "Neoceratodus_forsteri"
check <- name.check(full_data_tree, dat)
rm_phy <- check$tree_not_data
rm_dat <- check$data_not_tree
pruned_data_tree <- drop.tip(full_data_tree, rm_phy)
pruned_dat <- subset(dat, subset = dat$Scientific_name %in% pruned_data_tree$tip, select = names(dat))
str(pruned_dat)
> 'data.frame': 142 obs. of  18 variables:
>  $ Common.name               : chr  "Adelie penguin" "Alpine swift " "American bald eagle " "American flamengo" ...
>  $ Scientific_name           : chr  "Pygoscelis_adeliae" "Tachymarptis_melba" "Haliaeetus_leucocephalus" "Phoenicopterus_ruber" ...
>  $ Domain                    : chr  "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
>  $ Kingdom                   : chr  "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
>  $ Phylum                    : chr  "Chordata" "Chordata" "Chordata" "Chordata" ...
>  $ Class                     : chr  "Aves" "Aves" "Aves" "Aves" ...
>  $ Order                     : chr  "Sphenisciformes" "Apodiformes" "Accipitriformes" "Phoenicopteriformes" ...
>  $ Family                    : chr  "Spheniscidae" "Apodidae" "Accipitridae" "Phoenicopteridae" ...
>  $ Genus                     : chr  "Pygoscelis" "Tachymarptis" "Haliaeetus" "Phoenicopterus" ...
>  $ Species                   : chr  "adeliae" "melba" "leucocephalus" "ruber" ...
>  $ Endo_ectotherm            : chr  "endo" "endo" "endo" "endo" ...
>  $ Adult_mass_grams          : num  5020 103 3175 3066 117 ...
>  $ log_mass                  : num  8.52 4.64 8.06 8.03 4.77 ...
>  $ Lifespan_years            : num  20 26 48 44 17 25 16 28.2 29.9 28.5 ...
>  $ Loglifespan               : num  1.3 1.41 1.68 1.64 1.23 ...
>  $ Average_Telomere_Length_kb: num  7.35 18.78 2.1 18.31 30.2 ...
>  $ logTL                     : num  0.866 1.274 0.322 1.263 1.48 ...
>  $ Telomerase_activity       : int  NA NA NA NA NA NA NA NA NA NA ...
head(pruned_dat)
>                                   Common.name          Scientific_name
> Pygoscelis_adeliae             Adelie penguin       Pygoscelis_adeliae
> Tachymarptis_melba              Alpine swift        Tachymarptis_melba
> Haliaeetus_leucocephalus American bald eagle  Haliaeetus_leucocephalus
> Phoenicopterus_ruber        American flamengo     Phoenicopterus_ruber
> Falco_sparverius           American kestrel           Falco_sparverius
> Ichthyaetus_audouinii          Andouin's gull    Ichthyaetus_audouinii
>                             Domain Kingdom   Phylum Class               Order
> Pygoscelis_adeliae       Eukaryota Metazoa Chordata  Aves     Sphenisciformes
> Tachymarptis_melba       Eukaryota Metazoa Chordata  Aves         Apodiformes
> Haliaeetus_leucocephalus Eukaryota Metazoa Chordata  Aves     Accipitriformes
> Phoenicopterus_ruber     Eukaryota Metazoa Chordata  Aves Phoenicopteriformes
> Falco_sparverius         Eukaryota Metazoa Chordata  Aves       Falconiformes
> Ichthyaetus_audouinii    Eukaryota Metazoa Chordata  Aves     Charadriiformes
>                                    Family          Genus       Species
> Pygoscelis_adeliae           Spheniscidae     Pygoscelis       adeliae
> Tachymarptis_melba               Apodidae   Tachymarptis         melba
> Haliaeetus_leucocephalus     Accipitridae     Haliaeetus leucocephalus
> Phoenicopterus_ruber     Phoenicopteridae Phoenicopterus         ruber
> Falco_sparverius               Falconidae          Falco    sparverius
> Ichthyaetus_audouinii             Laridae    Ichthyaetus     audouinii
>                          Endo_ectotherm Adult_mass_grams log_mass
> Pygoscelis_adeliae                 endo           5020.0 8.521384
> Tachymarptis_melba                 endo            102.7 4.641502
> Haliaeetus_leucocephalus           endo           3175.0 8.063378
> Phoenicopterus_ruber               endo           3066.0 8.028455
> Falco_sparverius                   endo            117.0 4.770685
> Ichthyaetus_audouinii              endo            535.0 6.284134
>                          Lifespan_years Loglifespan Average_Telomere_Length_kb
> Pygoscelis_adeliae                   20    1.301030                      7.350
> Tachymarptis_melba                   26    1.414973                     18.783
> Haliaeetus_leucocephalus             48    1.681241                      2.100
> Phoenicopterus_ruber                 44    1.643453                     18.310
> Falco_sparverius                     17    1.230449                     30.200
> Ichthyaetus_audouinii                25    1.397940                     26.550
>                              logTL Telomerase_activity
> Pygoscelis_adeliae       0.8662873                  NA
> Tachymarptis_melba       1.2737650                  NA
> Haliaeetus_leucocephalus 0.3222193                  NA
> Phoenicopterus_ruber     1.2626883                  NA
> Falco_sparverius         1.4800069                  NA
> Ichthyaetus_audouinii    1.4240645                  NA
pruned_data_tree
> 
> Phylogenetic tree with 142 tips and 141 internal nodes.
> 
> Tip labels:
>   Torpedo_marmorata, Urolophus_paucimaculatus, Squalus_acanthias, Mustelus_asterias, Callorhinchus_milii, Oncorhynchus_mykiss, ...
> Node labels:
>   , '36', '37', '14', '25', '286', ...
> 
> Rooted; includes branch lengths.
name.check(pruned_data_tree, pruned_dat)
> [1] "OK"
hist(pruned_dat$Lifespan_years, main = "raw variable")

hist(log1p(pruned_dat$Lifespan_years), main = "log-transformed")

names(pruned_dat)
>  [1] "Common.name"                "Scientific_name"           
>  [3] "Domain"                     "Kingdom"                   
>  [5] "Phylum"                     "Class"                     
>  [7] "Order"                      "Family"                    
>  [9] "Genus"                      "Species"                   
> [11] "Endo_ectotherm"             "Adult_mass_grams"          
> [13] "log_mass"                   "Lifespan_years"            
> [15] "Loglifespan"                "Average_Telomere_Length_kb"
> [17] "logTL"                      "Telomerase_activity"
pruned_dat$log.lifespan <- log1p(pruned_dat$Lifespan_years)
pruned_dat$log.lifespan
>   [1] 3.0445224 3.2958369 3.8918203 3.8066625 2.8903718 3.2580965 2.8332133
>   [8] 3.3741687 3.4307562 3.3843903 3.1267605 3.4339872 3.5263605 3.3945084
>  [15] 3.7909847 2.7725887 2.6672282 3.7841896 2.7972813 3.7471484 4.1108739
>  [22] 3.6109179 3.5807373 2.4510051 2.0014800 3.9512437 3.5945688 3.1354942
>  [29] 3.4339872 3.4563167 2.4159138 3.1780538 2.4069451 2.9444390 3.8712010
>  [36] 3.2542430 3.4011974 3.2580965 2.5726122 3.9318256 2.3978953 2.6390573
>  [43] 2.5649494 3.2503745 2.8903718 3.9318256 1.7917595 2.6390573 2.3025851
>  [50] 3.7376696 3.8712010 1.8718022 1.7917595 1.7917595 0.2070142 0.9162907
>  [57] 1.6094379 3.2580965 2.7725887 2.0794415 2.3978953 3.0445224 2.9444390
>  [64] 2.4849066 3.1780538 2.9601051 1.9459101 2.5649494 4.3307333 3.0445224
>  [71] 4.4543473 2.7146947 5.3565863 4.3567088 3.7013020 3.1696856 3.1267605
>  [78] 3.0819100 3.9627161 3.3809947 3.3322045 3.1863526 2.9856819 3.0445224
>  [85] 3.6027768 3.7400477 3.3322045 3.8022081 3.4339872 2.9957323 3.3068867
>  [92] 2.9755296 3.5553481 2.5649494 3.3672958 2.8449094 2.0281482 2.5176965
>  [99] 1.4350845 2.7600099 2.5494452 2.1972246 2.6390573 2.0794415 2.2721259
> [106] 2.1860513 3.8286414 3.4657359 4.0604430 3.6243409 3.4657359 3.7376696
> [113] 4.8162412 3.8732822 3.4404181 4.0253517 4.0943446 4.2341065 4.1125119
> [120] 3.6454499 3.6532523 3.7135721 3.6888795 4.1896547 4.3894986 3.1945831
> [127] 2.7788193 3.2027464 1.9459101 3.4657359 1.5686159 1.6094379 2.5336968
> [134] 2.4932055 4.1896547 4.3567088 2.4849066 2.3978953 3.3250360 4.7957905
> [141] 3.4339872 4.1271344



Figure 1

#fonts()

#png("figure1.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure1.pdf")

par(family = "Monaco")
plot(NA, type = "n", xaxt = "n", ylab = "", xlab = "", xlim = c(0, 5), ylim = c(0, 5), axes = FALSE)
box()
text(x = 2.5, y = 5, "Body size", cex = 1.4)
text(x = 0.5, y = 0.5, "Lifespan", cex = 1.4)
text(x = 4.5, y = 0.5, "Endothermy", cex = 1.4)
text(x = 2.5, y = 2.5, "Telomere\nlength", cex = 1.4)

Arrows(2.1, 4.8, 0.7, 0.7, lwd = 2, code = 3, arr.adj = 2, arr.type="triangle", srt= 50, col = "gray")
Arrows(3, 4.8, 4.5, 0.7, lwd = 2, code = 3, arr.adj = 2, arr.type="triangle", srt= 50, col = "gray")
Arrows(1.2, 0.5, 3.5, 0.5, lwd = 2, code = 3, arr.adj = 2, arr.type="triangle", srt= 50, col = "gray")

Arrows(1, 0.7, 2.3, 2, lwd = 2, code = 2, arr.adj = 2, arr.type="triangle", srt= 50)
Arrows(3.9, 0.7, 2.7, 2, lwd = 2, code = 2, arr.adj = 2, arr.type="triangle", srt= 50)
Arrows(2.5, 4.8, 2.5, 2.8, lwd = 2, code = 2, arr.adj = 2, arr.type="triangle", srt= 50)

##dev.off()



Figure 2

#png("figure2.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure2.pdf")

plotTree(pruned_data_tree, fsize = 0.25, ftype = "i", mar = c(1, 1, 1, 3))

cladebox <- function(tree, node, spp, ysize, color,...){
    obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
    h <- max(nodeHeights(tree))
    parent <- tree$edge[which(tree$edge[ , 2] == node), 1]
    x0 <- max(c(obj$xx[node] + obj$xx[parent])/2, obj$xx[node]-0.01*h)
    x1 <- obj$x.lim[2]
    dd <- getDescendants(tree, node)
    y0 <- min(range(obj$yy[dd]))-0.6
    y1 <- max(range(obj$yy[dd]))+0.6
    cols <- colorRampPalette(c("white", color))(100)
    x0 <- seq(x0, x1, length.out = 101)[1:100]
    x1 <- seq(x0[1], x1, length.out = 101)[2:101]
    for(j in 1:100){
        polygon(c(x0[j], x1[j], x1[j], x0[j]), c(y0, y0, y1, y1), col = make.transparent(cols[j], 0.6),
                border = 0)
    }
    
    spp <- gsub("_", " ", as.character(spp))
    x2<- runif(100, obj$x.lim[2] + 10, obj$x.lim[2]+50)
    y3 <- seq(y0, y1, 5)
   
    for(i in spp){
        x2
        uuid <- get_uuid(name = i, n = 1)
        img <- get_phylopic(uuid = uuid)
        i <- gsub(" ", "_", i) ## This is actually not required as I think rphylopic can get spp names with "_"
        nodes <- sapply(i, grep, x = tree$tip.label)
        for(j in nodes){
            add_phylopic_base(img = img, x = sample(x2, 1), y = j, ysize = ysize, color = color, fill = color)
        }
       
    }

}

#nodelabels(cex = 0.5)

chel.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Reptilia"][8]
bird.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Aves"][c(4, 10, 11, 12, 27, 28, 5)]
alg.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Reptilia"][2]
squ.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Reptilia"][4]
mammal.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Mammalia"][c(4, 5, 2, 16, 30, 37, 47, 53, 55, 49)]
ost.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Fish"][c(7, 4, 8)]
chon.spp <- pruned_dat$Scientific_name[pruned_dat$Class == "Fish"][c(19, 24)]

cladebox(pruned_data_tree, 282, spp = chel.spp, 3, color = "#e93e3a")
cladebox(pruned_data_tree, 238, spp = bird.spp, 3, color = "#759c8a")
cladebox(pruned_data_tree, 237, spp = alg.spp, 2, color = "brown")
cladebox(pruned_data_tree, 233, spp = squ.spp, 1, color = "#e7dddf")
cladebox(pruned_data_tree, 170, spp = mammal.spp, 3, color = "#dd9f40")
cladebox(pruned_data_tree, 149, spp = ost.spp, 2, color = "#a39db8")
cladebox(pruned_data_tree, 144, spp = chon.spp, 2,  color = "#efe0a8")





#dev.off()

## PGLS models

names(pruned_dat)
>  [1] "Common.name"                "Scientific_name"           
>  [3] "Domain"                     "Kingdom"                   
>  [5] "Phylum"                     "Class"                     
>  [7] "Order"                      "Family"                    
>  [9] "Genus"                      "Species"                   
> [11] "Endo_ectotherm"             "Adult_mass_grams"          
> [13] "log_mass"                   "Lifespan_years"            
> [15] "Loglifespan"                "Average_Telomere_Length_kb"
> [17] "logTL"                      "Telomerase_activity"       
> [19] "log.lifespan"
pruned_dat$new.log.TL <- log(pruned_dat$Average_Telomere_Length_kb)
pruned_dat$new.log.TL
>   [1] 1.9947003 2.9329522 0.7419373 2.9074474 3.4078419 3.2790297 2.2586332
>   [8] 2.4638532 2.6803364 2.2828925 2.4024304 2.9806186 2.5494452 2.6253930
>  [15] 2.2813615 1.7101878 1.1631508 1.9586853 3.8971124 2.8534772 2.7507904
>  [22] 2.9460167 2.5257286 1.8373700 3.1527360 2.6567569 2.1488510 0.2623643
>  [29] 2.4680995 2.5877640 2.0918641 2.6946272 3.0568274 2.4535880 2.0893919
>  [36] 3.7471484 2.0412203 2.6100698 2.6100698 2.2925348 2.7948393 2.6318888
>  [43] 2.6173958 2.1587147 1.3737156 2.5257286 2.0149030 1.4011830 1.9459101
>  [50] 1.5260563 1.5303947 2.7725887 1.4469190 1.4469190 1.8885837 2.0502702
>  [57] 1.7917595 2.4313300 1.2354715 1.7155981 1.7011051 1.8148247 1.6601310
>  [64] 2.9957323 2.5006159 1.0986123 3.2023399 3.1941732 2.5257286 1.0986123
>  [71] 3.1871787 3.9120230 2.1972246 2.5649494 2.3025851 2.8903718 2.9856819
>  [78] 2.3418058 2.8332133 2.5649494 2.7080502 2.6390573 2.6390573 2.8903718
>  [85] 1.6094379 2.0794415 2.7080502 2.4849066 3.2695689 3.2188758 3.9120230
>  [92] 2.2874715 3.4011974 3.2580965 2.1972246 2.5649494 3.5553481 3.6375862
>  [99] 2.5329028 2.7080502 3.9120230 3.9120230 4.3820266 2.9957323 3.6109179
> [106] 3.6109179 2.3025851 2.4849066 2.9285235 2.4849066 2.4849066 2.8332133
> [113] 2.1972246 1.9459101 2.1972246 2.3025851 2.3025851 2.4510051 2.7343675
> [120] 2.9444390 2.7013612 2.7725887 2.7146947 2.6390573 2.7080502 1.9459101
> [127] 2.3025851 3.9120230 2.1972246 2.7725887 3.6888795 3.6888795 2.7080502
> [134] 3.1354942 2.7825391 3.4307562 2.9606231 2.9231616 3.2976870 2.9957323
> [141] 3.9120230 4.0943446
model1 <- gls(new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan + log_mass:Endo_ectotherm + log.lifespan:Endo_ectotherm + log_mass:log.lifespan:Endo_ectotherm, correlation = corBrownian(phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(model1)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan +      log_mass:Endo_ectotherm + log.lifespan:Endo_ectotherm + log_mass:log.lifespan:Endo_ectotherm 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   343.6862 370.2886 -162.8431
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                                              Value Std.Error    t-value p-value
> (Intercept)                               3.233397 1.3121680  2.4641640  0.0150
> log_mass                                 -0.034819 0.2118296 -0.1643744  0.8697
> log.lifespan                             -0.138426 0.3991831 -0.3467736  0.7293
> Endo_ectothermendo                       -0.175596 1.7674781 -0.0993482  0.9210
> log_mass:log.lifespan                     0.006902 0.0705437  0.0978350  0.9222
> log_mass:Endo_ectothermendo               0.011989 0.2326388  0.0515342  0.9590
> log.lifespan:Endo_ectothermendo           0.013390 0.4800741  0.0278912  0.9778
> log_mass:log.lifespan:Endo_ectothermendo -0.012800 0.0752810 -0.1700348  0.8652
> 
>  Correlation: 
>                                          (Intr) lg_mss lg.lfs End_ct lg_m:.
> log_mass                                 -0.469                            
> log.lifespan                             -0.615  0.535                     
> Endo_ectothermendo                       -0.495  0.279  0.453              
> log_mass:log.lifespan                     0.456 -0.925 -0.730 -0.279       
> log_mass:Endo_ectothermendo               0.427 -0.911 -0.487 -0.418  0.842
> log.lifespan:Endo_ectothermendo           0.512 -0.445 -0.832 -0.623  0.607
> log_mass:log.lifespan:Endo_ectothermendo -0.427  0.867  0.684  0.405 -0.937
>                                          lg_:E_ lg.:E_
> log_mass                                              
> log.lifespan                                          
> Endo_ectothermendo                                    
> log_mass:log.lifespan                                 
> log_mass:Endo_ectothermendo                           
> log.lifespan:Endo_ectothermendo           0.562       
> log_mass:log.lifespan:Endo_ectothermendo -0.922 -0.732
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.00833399 -0.09146447  0.13823694  0.30407049  0.92539247 
> 
> Residual standard error: 2.098608 
> Degrees of freedom: 142 total; 134 residual
model2 <- update(model1, ~ .-log_mass:log.lifespan:Endo_ectotherm)
summary(model2)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan +      log_mass:Endo_ectotherm + log.lifespan:Endo_ectotherm 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   341.7168 365.3634 -162.8584
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                                      Value Std.Error    t-value p-value
> (Intercept)                      3.1380986 1.1821750  2.6545128  0.0089
> log_mass                        -0.0035886 0.1051473 -0.0341291  0.9728
> log.lifespan                    -0.0920198 0.2902581 -0.3170273  0.7517
> Endo_ectothermendo              -0.0538850 1.6102252 -0.0334643  0.9734
> log_mass:log.lifespan           -0.0043401 0.0245142 -0.1770442  0.8597
> log_mass:Endo_ectothermendo     -0.0245012 0.0894885 -0.2737915  0.7847
> log.lifespan:Endo_ectothermendo -0.0463691 0.3258585 -0.1422984  0.8871
> 
>  Correlation: 
>                                 (Intr) lg_mss lg.lfs End_ct lg_m:. lg_:E_
> log_mass                        -0.219                                   
> log.lifespan                    -0.489 -0.160                            
> Endo_ectothermendo              -0.389 -0.159  0.264                     
> log_mass:log.lifespan            0.175 -0.648 -0.349  0.317              
> log_mass:Endo_ectothermendo      0.095 -0.575  0.510 -0.125 -0.165       
> log.lifespan:Endo_ectothermendo  0.324  0.558 -0.666 -0.525 -0.332 -0.430
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -0.99962915 -0.09501665  0.14604616  0.32250953  0.93359053 
> 
> Residual standard error: 2.098834 
> Degrees of freedom: 142 total; 135 residual
model3 <- update(model2, ~ .-log.lifespan:Endo_ectotherm)
summary(model3)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan +      log_mass:Endo_ectotherm 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   339.7381 360.4289 -162.8691
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                                 Value Std.Error    t-value p-value
> (Intercept)                  3.192564 1.1144614  2.8646697  0.0048
> log_mass                     0.004768 0.0869066  0.0548589  0.9563
> log.lifespan                -0.119529 0.2157295 -0.5540687  0.5804
> Endo_ectothermendo          -0.174090 1.3659091 -0.1274533  0.8988
> log_mass:log.lifespan       -0.005499 0.0230394 -0.2386600  0.8117
> log_mass:Endo_ectothermendo -0.029972 0.0805158 -0.3722555  0.7103
> 
>  Correlation: 
>                             (Intr) lg_mss lg.lfs End_ct lg_m:.
> log_mass                    -0.509                            
> log.lifespan                -0.387  0.343                     
> Endo_ectothermendo          -0.272  0.190 -0.135              
> log_mass:log.lifespan        0.317 -0.591 -0.810  0.178       
> log_mass:Endo_ectothermendo  0.273 -0.448  0.333 -0.456 -0.361
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -0.99393208 -0.09318428  0.14374430  0.32634461  0.94174318 
> 
> Residual standard error: 2.098991 
> Degrees of freedom: 142 total; 136 residual
model4 <- update(model3, ~ .-log_mass:log.lifespan)
summary(model4)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:Endo_ectotherm 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   337.7976 355.5325 -162.8988
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                                 Value Std.Error    t-value p-value
> (Intercept)                  3.276795 1.0534568  3.1105164  0.0023
> log_mass                    -0.007483 0.0698850 -0.1070806  0.9149
> log.lifespan                -0.161249 0.1259793 -1.2799652  0.2027
> Endo_ectothermendo          -0.116185 1.3395539 -0.0867343  0.9310
> log_mass:Endo_ectothermendo -0.036908 0.0748303 -0.4932183  0.6226
> 
>  Correlation: 
>                             (Intr) lg_mss lg.lfs End_ct
> log_mass                    -0.421                     
> log.lifespan                -0.235 -0.286              
> Endo_ectothermendo          -0.352  0.371  0.015       
> log_mass:Endo_ectothermendo  0.438 -0.878  0.074 -0.427
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -0.99267537 -0.09429773  0.14396653  0.32720827  0.94298031 
> 
> Residual standard error: 2.099431 
> Degrees of freedom: 142 total; 137 residual
model5 <- update(model4, ~ .-log_mass:Endo_ectotherm)
summary(model5)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   336.0495 350.8286 -163.0247
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                        Value Std.Error   t-value p-value
> (Intercept)         3.504551 0.9442549  3.711446  0.0003
> log_mass           -0.037761 0.0333064 -1.133747  0.2589
> log.lifespan       -0.156676 0.1252926 -1.250477  0.2132
> Endo_ectothermendo -0.398370 1.2079020 -0.329803  0.7420
> 
>  Correlation: 
>                    (Intr) lg_mss lg.lfs
> log_mass           -0.084              
> log.lifespan       -0.298 -0.465       
> Endo_ectothermendo -0.203 -0.009  0.052
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -0.99465171 -0.09075665  0.13924989  0.31601039  0.93865079 
> 
> Residual standard error: 2.101294 
> Degrees of freedom: 142 total; 138 residual
model6 <- update(model5, ~ .-Endo_ectotherm)
summary(model6)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   334.1614 345.9847 -163.0807
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                  Value Std.Error   t-value p-value
> (Intercept)   3.441364 0.9216444  3.733939  0.0003
> log_mass     -0.037861 0.0331981 -1.140468  0.2561
> log.lifespan -0.154544 0.1247240 -1.239089  0.2174
> 
>  Correlation: 
>              (Intr) lg_mss
> log_mass     -0.088       
> log.lifespan -0.294 -0.466
> 
> Standardized residuals:
>           Min            Q1           Med            Q3           Max 
> -1.1565570775 -0.2145557854  0.0002704793  0.1739263041  0.7765132875 
> 
> Residual standard error: 2.102122 
> Degrees of freedom: 142 total; 139 residual
model7 <- update(model6, ~ .-log_mass)
summary(model7)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log.lifespan 
>   Data: pruned_dat 
>        AIC      BIC   logLik
>   333.4839 342.3514 -163.742
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                  Value Std.Error   t-value p-value
> (Intercept)   3.348911 0.9190577  3.643853  0.0004
> log.lifespan -0.220773 0.1104984 -1.997977  0.0477
> 
>  Correlation: 
>              (Intr)
> log.lifespan -0.38 
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.13370657 -0.20204102 -0.01077654  0.15669655  0.78439669 
> 
> Residual standard error: 2.111934 
> Degrees of freedom: 142 total; 140 residual
## Model diagnosis


layout(matrix(c(0, 0, 0, 0,
                1, 1, 2, 2,
                1, 1, 2, 2,
                0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))



## Checking homogeneity of variance


plot(fitted(model7), resid(model7), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals")
abline(h = 0, col = "darkorange", lwd = 2)


## Checking normality

qqnorm(resid(model7), col = "darkgrey")
qqline(resid(model7), col = "dodgerblue", lwd = 2)



Figure 5

#png("figure3.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure3.pdf")


par(mar = c(5, 4, 1, 4))

plot(new.log.TL ~ log.lifespan , data = pruned_dat, pch = 21, bg = c(alpha("purple", 0.5), alpha("orange", 0.5), alpha("pink", 0.7), alpha("gold", 0.5))[as.numeric(as.factor(pruned_dat$Class))], col = c(alpha("purple", 0.5), alpha("orange", 0.5), alpha("pink", 0.7), alpha("gold", 0.5))[as.numeric(as.factor(pruned_dat$Class))], las = 1, xlab = "Log lifespan (yrs)", ylab = "Log telomere length (kb)", type = "n")

grid(nx = NULL, ny = NULL, col = "lightgray", lwd = 1)

par(new = TRUE)

plot(new.log.TL ~ log.lifespan , data = pruned_dat, pch = 21, bg = c(alpha("purple", 0.5), alpha("orange", 0.5), alpha("pink", 0.7), alpha("gold", 0.5))[as.numeric(as.factor(pruned_dat$Class))], col = c(alpha("purple", 0.5), alpha("orange", 0.5), alpha("pink", 0.7), alpha("gold", 0.5))[as.numeric(as.factor(pruned_dat$Class))], las = 1, xlab = "Log lifespan (yrs)", ylab = "Log telomere length (kb)")

par(xpd = TRUE)

uuid.aves <- get_uuid(name = "Accipiter_gentilis", n = 1)
img.aves <- get_phylopic(uuid = uuid.aves)
add_phylopic_base(img = img.aves, x = 6, y = 2.4, ysize = 0.1, col = alpha("purple", 0.5), fill = alpha("purple", 0.5))

uuid.fish <- get_uuid(name = "Oreochromis niloticus", n = 1)
img.fish <- get_phylopic(uuid = uuid.fish)
add_phylopic_base(img = img.fish, x = 6, y = 2.2, ysize = 0.09, col = alpha("orange", 0.5), fill = alpha("orange", 0.5))

uuid.mammals <- get_uuid(name = "Elephas_maximus", n = 1)
img.mammals <- get_phylopic(uuid = uuid.mammals)
add_phylopic_base(img = img.mammals, x = 6, y = 2, ysize = 0.2, col = alpha("pink", 0.7), fill = alpha("pink", 0.7), horizontal = TRUE)

uuid.reptilia <- get_uuid(name = "Alligator_mississippiensis", n = 1)
img.reptilia <- get_phylopic(uuid = uuid.reptilia)
add_phylopic_base(img = img.reptilia, x = 6, y = 1.8, ysize = 0.1, col = alpha("gold", 0.5), fill = alpha("gold", 0.5))


SSX <- sum(round((pruned_dat$new.log.TL - mean(pruned_dat$new.log.TL))^2), 2)
s2 <- var(pruned_dat$new.log.TL)
n <- length(pruned_dat$new.log.TL)
x <- seq(min(pruned_dat$log.lifespan), max(pruned_dat$log.lifespan), 0.06)
m.x <- mean(round(pruned_dat$new.log.TL, 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(model7)[1] + coef(model7)[2]*x) + ic.s
lower.i <- (coef(model7)[1] + coef(model7)[2]*x) + ic.i

##par(new = TRUE)

lines(x = x, y = (coef(model7)[1] + (coef(model7)[2]*x)), lwd = 2, col = alpha("black", 0.5))

polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("gold", 0.3))

Figure 4

## Phylogenetic signal (Blomberg et al.'s K statitics

telo.length <- setNames(pruned_dat$new.log.TL, rownames(pruned_dat))
telo.length
>         Pygoscelis_adeliae         Tachymarptis_melba 
>                  1.9947003                  2.9329522 
>   Haliaeetus_leucocephalus       Phoenicopterus_ruber 
>                  0.7419373                  2.9074474 
>           Falco_sparverius      Ichthyaetus_audouinii 
>                  3.4078419                  3.2790297 
>            Hirundo_rustica           Branta_leucopsis 
>                  2.2586332                  2.4638532 
>             Cepphus_grylle           Rissa_tridactyla 
>                  2.6803364                  2.2828925 
>              Turdus_merula              Gallus_gallus 
>                  2.4024304                  2.9806186 
>             Sterna_hirundo            Calidris_alpina 
>                  2.5494452                  2.6253930 
>      Haematopus_ostralegus    Aphelocoma_coerulescens 
>                  2.2813615                  1.7101878 
>        Chrysolophus_pictus              Fregata_minor 
>                  1.1631508                  1.9586853 
>                Parus_major                Gyps_fulvus 
>                  3.8971124                  2.8534772 
>        Strigops_habroptila      Oceanodroma_leucorhoa 
>                  2.7507904                  2.9460167 
>               Larus_fuscus      Tachycineta_albilinea 
>                  2.5257286                  1.8373700 
>        Colinus_virginianus         Fulmarus_glacialis 
>                  3.1527360                  2.6567569 
>          Macronectes_halli         Accipiter_gentilis 
>                  2.1488510                  0.2623643 
>          Amazona_amazonica          Buteo_jamaicensis 
>                  2.4680995                  2.5877640 
>            Vireo_olivaceus                  Sula_sula 
>                  2.0918641                  2.6946272 
>            Riparia_riparia  Passerculus_sandwichensis 
>                  3.0568274                  2.4535880 
>      Macronectes_giganteus            Buteo_swainsoni 
>                  2.0893919                  3.7471484 
>                Uria_lomvia     Aphelocoma_ultramarina 
>                  2.0412203                  2.6100698 
>        Tachycineta_bicolor           Diomedea_exulans 
>                  2.6100698                  2.2925348 
>           Lonchura_striata        Meleagris_gallopavo 
>                  2.7948393                  2.6318888 
>        Taeniopygia_guttata       Calonectris_diomedea 
>                  2.6173958                  2.1587147 
>  Acrocephalus_sechellensis          Anguilla_rostrata 
>                  1.3737156                  2.5257286 
>            Oryzias_latipes       Pseudocaranx_wrighti 
>                  2.0149030                  1.4011830 
>      Oreochromis_niloticus          Carassius_auratus 
>                  1.9459101                  1.5260563 
>            Cyprinus_carpio                Danio_rerio 
>                  1.5303947                  2.7725887 
>      Xiphophorus_maculatus       Xiphophorus_hellerii 
>                  1.4469190                  1.4469190 
>     Nothobranchius_furzeri    Nothobranchius_rachovii 
>                  1.8885837                  2.0502702 
>      Fundulus_heteroclitus               Gadus_morhua 
>                  1.7917595                  2.4313300 
>       Dicentrarchus_labrax   Acanthopagrus_schlegelii 
>                  1.2354715                  1.7155981 
>     Upeneichthys_vlamingii          Macquaria_ambigua 
>                  1.7011051                  1.8148247 
>  Lutjanus_argentimaculatus        Oncorhynchus_mykiss 
>                  1.6601310                  2.9957323 
>    Platycephalus_bassensis          Mustelus_asterias 
>                  2.5006159                  1.0986123 
>        Callorhinchus_milii   Urolophus_paucimaculatus 
>                  3.2023399                  3.1941732 
>          Squalus_acanthias          Torpedo_marmorata 
>                  2.5257286                  1.0986123 
>      Neoceratodus_forsteri            Setifer_setosus 
>                  3.1871787                  3.9120230 
>         Balaena_mysticetus      Eschrichtius_robustus 
>                  2.1972246                  2.5649494 
>     Giraffa_camelopardalis                 Ovis_aries 
>                  2.3025851                  2.8903718 
>          Rangifer_tarandus               Capra_hircus 
>                  2.9856819                  2.3418058 
>         Tursiops_truncatus        Camelus_dromedarius 
>                  2.8332133                  2.5649494 
>                 Sus_scrofa          Muntiacus_reevesi 
>                  2.7080502                  2.6390573 
>          Muntiacus_muntjak                 Bos_taurus 
>                  2.6390573                  2.8903718 
>     Zalophus_californianus            Crocuta_crocuta 
>                  1.6094379                  2.0794415 
>                Canis_lupus            Ursus_maritimus 
>                  2.7080502                  2.4849066 
>                Felis_catus            Ailurus_fulgens 
>                  3.2695689                  3.2188758 
>            Panthera_tigris                Meles_meles 
>                  3.9120230                  2.2874715 
>           Myotis_lucifugus      Tadarida_brasiliensis 
>                  3.4011974                  3.2580965 
>       Pteropus_rodricensis  Chaetophractus_vellerosus 
>                  2.1972246                  2.5649494 
>       Didelphis_virginiana       Atelerix_albiventris 
>                  3.5553481                  3.6375862 
>              Sorex_araneus          Procavia_capensis 
>                  2.5329028                  2.7080502 
>         Lepus_californicus       Sylvilagus_aquaticus 
>                  3.9120230                  3.9120230 
>      Oryctolagus_cuniculus          Ochotona_princeps 
>                  4.3820266                  2.9957323 
> Macroscelides_proboscideus       Galegeeska_rufescens 
>                  3.6109179                  3.6109179 
>        Ceratotherium_simum               Equus_grevyi 
>                  2.3025851                  2.4849066 
>             Equus_caballus            Tapirus_indicus 
>                  2.9285235                  2.4849066 
>    Myrmecophaga_tridactyla        Choloepus_hoffmanni 
>                  2.4849066                  2.8332133 
>               Homo_sapiens           Ateles_geoffroyi 
>                  2.1972246                  1.9459101 
>           Saimiri_sciureus               Pan_paniscus 
>                  2.1972246                  2.3025851 
>             Pongo_pygmaeus            Pan_troglodytes 
>                  2.3025851                  2.4510051 
>            Gorilla_gorilla                Lemur_catta 
>                  2.7343675                  2.9444390 
>          Macaca_nemestrina             Macaca_mulatta 
>                  2.7013612                  2.7725887 
>        Macaca_fascicularis         Loxodonta_africana 
>                  2.7146947                  2.6390573 
>            Elephas_maximus          Castor_canadensis 
>                  2.7080502                  1.9459101 
>  Hydrochoerus_hydrochaeris       Sciurus_carolinensis 
>                  2.3025851                  3.9120230 
>            Aplodontia_rufa      Heterocephalus_glaber 
>                  2.1972246                  2.7725887 
>          Rattus_norvegicus               Mus_musculus 
>                  3.6888795                  3.6888795 
>                Tupaia_tana           Tupaia_belangeri 
>                  2.7080502                  3.1354942 
>         Alligator_sinensis Alligator_mississippiensis 
>                  2.7825391                  3.4307562 
>           Zootoca_vivipara             Lacerta_agilis 
>                  2.9606231                  2.9231616 
>              Liasis_fuscus           Emys_orbicularis 
>                  3.2976870                  2.9957323 
>          Trachemys_scripta            Chrysemys_picta 
>                  3.9120230                  4.0943446
k_tl <- phylosig(pruned_data_tree, telo.length, test = TRUE, nsim = 10000)
attributes(k_tl)
> $names
> [1] "K"     "P"     "sim.K"
> 
> $class
> [1] "phylosig"
> 
> $method
> [1] "K"
> 
> $test
> [1] TRUE
> 
> $se
> [1] FALSE
head(k_tl$sim.K)
> [1] 0.13086102 0.05921165 0.05554814 0.03498161 0.03009382 0.03544119
k_tl$K
> [1] 0.130861
k_tl$P
> [1] 2e-04
#png("figure4.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure4.pdf")

hist(k_tl$sim.K, bty = "o", ylim = c(0, 3000), las = 1, ylab = "Null distribution of K", xlab = "K", main = "")
abline(v = k_tl$K, lwd = 2, lty = "dotted")
text(x = k_tl$K - 0.02, y = 2500, "Observed value \n of K")
box()

#dev.off()



## PGLS models within classes

## Mammals

pruned_tree_mammals <- drop.tip(pruned_data_tree, pruned_dat$Scientific_name[!pruned_dat$Class == "Mammalia"])

name.check(pruned_tree_mammals, data = pruned_dat[pruned_dat$Class == "Mammalia", ])
> [1] "OK"
model8 <- gls(new.log.TL ~ log_mass + log.lifespan + log_mass:log.lifespan, correlation = corBrownian(phy = pruned_tree_mammals, form = ~Scientific_name), data = pruned_dat[pruned_dat$Class == "Mammalia", ], method = "ML")
summary(model8)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan + log_mass:log.lifespan 
>   Data: pruned_dat[pruned_dat$Class == "Mammalia", ] 
>        AIC      BIC    logLik
>   96.58006 107.2957 -43.29003
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                           Value Std.Error   t-value p-value
> (Intercept)            3.316300 0.8970282  3.696985  0.0005
> log_mass               0.001644 0.0777088  0.021161  0.9832
> log.lifespan           0.003464 0.2564737  0.013504  0.9893
> log_mass:log.lifespan -0.010498 0.0215340 -0.487518  0.6277
> 
>  Correlation: 
>                       (Intr) lg_mss lg.lfs
> log_mass              -0.712              
> log.lifespan          -0.778  0.682       
> log_mass:log.lifespan  0.745 -0.920 -0.853
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.45577217 -0.57238601 -0.26650191  0.08909956  1.42213153 
> 
> Residual standard error: 0.8803278 
> Degrees of freedom: 63 total; 59 residual
model9 <- update(model8, ~ .-log_mass:log.lifespan)
summary(model9)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan 
>   Data: pruned_dat[pruned_dat$Class == "Mammalia", ] 
>        AIC      BIC    logLik
>   94.83334 103.4059 -43.41667
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                  Value Std.Error   t-value p-value
> (Intercept)   3.642223 0.5942858  6.128739  0.0000
> log_mass     -0.033205 0.0302842 -1.096433  0.2773
> log.lifespan -0.103183 0.1330316 -0.775630  0.4410
> 
>  Correlation: 
>              (Intr) lg_mss
> log_mass     -0.102       
> log.lifespan -0.410 -0.500
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.42754826 -0.56875973 -0.29030500  0.09795673  1.42956182 
> 
> Residual standard error: 0.8820992 
> Degrees of freedom: 63 total; 60 residual
model10 <- update(model9, ~ .-log.lifespan)
summary(model10)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass 
>   Data: pruned_dat[pruned_dat$Class == "Mammalia", ] 
>        AIC      BIC    logLik
>   93.46188 99.89128 -43.73094
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                 Value Std.Error  t-value p-value
> (Intercept)  3.453449 0.5403901  6.39066  0.0000
> log_mass    -0.044940 0.0261481 -1.71867  0.0907
> 
>  Correlation: 
>          (Intr)
> log_mass -0.388
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -1.4666567 -0.5984019 -0.3109612  0.0536049  1.4274522 
> 
> Residual standard error: 0.8865104 
> Degrees of freedom: 63 total; 61 residual
## Birds + Reptiles

bird.reptiles <- rbind(pruned_dat[pruned_dat$Class == "Aves", ], pruned_dat[pruned_dat$Class == "Reptilia", ])

pruned_dat_minus_bird.reptiles <- pruned_dat[!pruned_dat$Scientific_name %in% bird.reptiles$Scientific_name, ]

pruned_tree_bird.reptiles <- drop.tip(pruned_data_tree, pruned_dat_minus_bird.reptiles$Scientific_name)
name.check(pruned_tree_bird.reptiles, bird.reptiles)
> [1] "OK"
model11 <- gls(new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan + log_mass:Endo_ectotherm + log.lifespan:Endo_ectotherm + log_mass:log.lifespan:Endo_ectotherm, correlation = corBrownian(phy = pruned_tree_bird.reptiles, form = ~Scientific_name), data = bird.reptiles, method = "ML")
summary(model11)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan +      log_mass:Endo_ectotherm + log.lifespan:Endo_ectotherm + log_mass:log.lifespan:Endo_ectotherm 
>   Data: bird.reptiles 
>        AIC     BIC    logLik
>   162.1094 179.842 -72.05468
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                                              Value Std.Error    t-value p-value
> (Intercept)                               4.881621  7.254073  0.6729490  0.5044
> log_mass                                 -0.118243  1.236757 -0.0956076  0.9243
> log.lifespan                             -0.673792  2.146858 -0.3138505  0.7551
> Endo_ectothermendo                       -3.252267  8.043313 -0.4043442  0.6879
> log_mass:log.lifespan                     0.059966  0.347826  0.1724018  0.8639
> log_mass:Endo_ectothermendo               0.410976  1.317415  0.3119564  0.7565
> log.lifespan:Endo_ectothermendo           1.217913  2.343425  0.5197150  0.6058
> log_mass:log.lifespan:Endo_ectothermendo -0.190072  0.375291 -0.5064671  0.6150
> 
>  Correlation: 
>                                          (Intr) lg_mss lg.lfs End_ct lg_m:.
> log_mass                                 -0.900                            
> log.lifespan                             -0.951  0.909                     
> Endo_ectothermendo                       -0.897  0.807  0.854              
> log_mass:log.lifespan                     0.871 -0.979 -0.931 -0.778       
> log_mass:Endo_ectothermendo               0.845 -0.939 -0.853 -0.873  0.920
> log.lifespan:Endo_ectothermendo           0.871 -0.833 -0.916 -0.922  0.853
> log_mass:log.lifespan:Endo_ectothermendo -0.807  0.908  0.863  0.847 -0.927
>                                          lg_:E_ lg.:E_
> log_mass                                              
> log.lifespan                                          
> Endo_ectothermendo                                    
> log_mass:log.lifespan                                 
> log_mass:Endo_ectothermendo                           
> log.lifespan:Endo_ectothermendo           0.903       
> log_mass:log.lifespan:Endo_ectothermendo -0.977 -0.932
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.01611101 -0.18310917 -0.04275159  0.11800721  0.55779323 
> 
> Residual standard error: 2.236097 
> Degrees of freedom: 53 total; 45 residual
model12 <- update(model11, ~ .-log_mass:log.lifespan:Endo_ectotherm)
summary(model12)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan +      log_mass:Endo_ectotherm + log.lifespan:Endo_ectotherm 
>   Data: bird.reptiles 
>        AIC      BIC   logLik
>   160.4106 176.1729 -72.2053
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                                      Value Std.Error    t-value p-value
> (Intercept)                      1.9158767  4.246819  0.4511321  0.6540
> log_mass                         0.4503777  0.514506  0.8753592  0.3859
> log.lifespan                     0.2646984  1.075332  0.2461552  0.8067
> Endo_ectothermendo               0.1969454  4.244780  0.0463971  0.9632
> log_mass:log.lifespan           -0.1033042  0.129554 -0.7973846  0.4293
> log_mass:Endo_ectothermendo     -0.2409376  0.278353 -0.8655822  0.3912
> log.lifespan:Endo_ectothermendo  0.1118884  0.843186  0.1326972  0.8950
> 
>  Correlation: 
>                                 (Intr) lg_mss lg.lfs End_ct lg_m:. lg_:E_
> log_mass                        -0.676                                   
> log.lifespan                    -0.852  0.593                            
> Endo_ectothermendo              -0.679  0.171  0.459                     
> log_mass:log.lifespan            0.554 -0.877 -0.693  0.034              
> log_mass:Endo_ectothermendo      0.449 -0.580 -0.094 -0.406  0.175       
> log.lifespan:Endo_ectothermendo  0.554  0.086 -0.610 -0.691 -0.077 -0.103
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -0.99688508 -0.16468067  0.01665759  0.12132122  0.57166703 
> 
> Residual standard error: 2.242461 
> Degrees of freedom: 53 total; 46 residual
model13 <- update(model12, ~ .-log.lifespan:Endo_ectotherm)
summary(model13)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:log.lifespan +      log_mass:Endo_ectotherm 
>   Data: bird.reptiles 
>        AIC      BIC    logLik
>   158.4309 172.2229 -72.21545
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                                  Value Std.Error    t-value p-value
> (Intercept)                  1.6035135  3.497602  0.4584608  0.6487
> log_mass                     0.4444743  0.507194  0.8763399  0.3853
> log.lifespan                 0.3517726  0.842963  0.4173048  0.6784
> Endo_ectothermendo           0.5859354  3.037761  0.1928840  0.8479
> log_mass:log.lifespan       -0.1019766  0.127810 -0.7978772  0.4290
> log_mass:Endo_ectothermendo -0.2371302  0.273962 -0.8655599  0.3911
> 
>  Correlation: 
>                             (Intr) lg_mss lg.lfs End_ct lg_m:.
> log_mass                    -0.874                            
> log.lifespan                -0.779  0.818                     
> Endo_ectothermendo          -0.492  0.320  0.065              
> log_mass:log.lifespan        0.719 -0.876 -0.936 -0.026       
> log_mass:Endo_ectothermendo  0.611 -0.576 -0.199 -0.664  0.168
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.00257514 -0.17211069 -0.01640975  0.10795757  0.56656400 
> 
> Residual standard error: 2.24289 
> Degrees of freedom: 53 total; 47 residual
model14 <- update(model13, ~ .-log_mass:log.lifespan)
summary(model14)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm + log_mass:Endo_ectotherm 
>   Data: bird.reptiles 
>       AIC      BIC    logLik
>   157.144 168.9657 -72.57198
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                                 Value Std.Error    t-value p-value
> (Intercept)                  3.610949 2.4204177  1.4918700  0.1423
> log_mass                     0.089910 0.2435472  0.3691674  0.7136
> log.lifespan                -0.278043 0.2946686 -0.9435777  0.3501
> Endo_ectothermendo           0.522535 3.0252042  0.1727271  0.8636
> log_mass:Endo_ectothermendo -0.200409 0.2690439 -0.7448938  0.4600
> 
>  Correlation: 
>                             (Intr) lg_mss lg.lfs End_ct
> log_mass                    -0.727                     
> log.lifespan                -0.433 -0.016              
> Endo_ectothermendo          -0.682  0.617  0.116       
> log_mass:Endo_ectothermendo  0.716 -0.903 -0.120 -0.669
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -0.98809648 -0.17093673 -0.03260668  0.10129896  0.56675072 
> 
> Residual standard error: 2.258029 
> Degrees of freedom: 53 total; 48 residual
model15 <- update(model14, ~ .-log_mass:Endo_ectotherm)
summary(model15)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan + Endo_ectotherm 
>   Data: bird.reptiles 
>        AIC      BIC    logLik
>   155.7531 165.6046 -72.87655
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                        Value Std.Error    t-value p-value
> (Intercept)         4.902176 1.6815740  2.9152307  0.0053
> log_mass           -0.073915 0.1041454 -0.7097328  0.4812
> log.lifespan       -0.304342 0.2912141 -1.0450783  0.3011
> Endo_ectothermendo -0.985788 2.2373682 -0.4406018  0.6614
> 
>  Correlation: 
>                    (Intr) lg_mss lg.lfs
> log_mass           -0.267              
> log.lifespan       -0.500 -0.292       
> Endo_ectothermendo -0.390  0.038  0.049
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -0.96251802 -0.14477137  0.01596402  0.15970774  0.58580223 
> 
> Residual standard error: 2.271042 
> Degrees of freedom: 53 total; 49 residual
model16 <- update(model15, ~ .-Endo_ectotherm)
summary(model16)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan 
>   Data: bird.reptiles 
>        AIC      BIC    logLik
>   153.9627 161.8438 -72.98133
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                  Value Std.Error    t-value p-value
> (Intercept)   4.613051 1.5357249  3.0038264  0.0042
> log_mass     -0.072153 0.1032265 -0.6989782  0.4878
> log.lifespan -0.298088 0.2885145 -1.0331819  0.3065
> 
>  Correlation: 
>              (Intr) lg_mss
> log_mass     -0.274       
> log.lifespan -0.523 -0.294
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.28077125 -0.44823059 -0.29285806 -0.05635584  0.50043747 
> 
> Residual standard error: 2.275537 
> Degrees of freedom: 53 total; 50 residual
model17 <- update(model16, ~ .-log_mass)
summary(model17)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log.lifespan 
>   Data: bird.reptiles 
>       AIC      BIC    logLik
>   152.478 158.3889 -73.23902
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                  Value Std.Error   t-value p-value
> (Intercept)   4.318816 1.4694824  2.939005  0.0049
> log.lifespan -0.357471 0.2743365 -1.303038  0.1984
> 
>  Correlation: 
>              (Intr)
> log.lifespan -0.657
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.28381370 -0.36979441 -0.23548939 -0.04838853  0.54703250 
> 
> Residual standard error: 2.286627 
> Degrees of freedom: 53 total; 51 residual
## Fish

pruned_tree_fish <- drop.tip(pruned_data_tree, pruned_dat$Scientific_name[!pruned_dat$Class == "Fish"])
name.check(pruned_tree_fish, data = pruned_dat[pruned_dat$Class == "Fish", ])
> [1] "OK"
model18 <- gls(new.log.TL ~ log_mass + log.lifespan + log_mass:log.lifespan, correlation = corBrownian(phy = pruned_tree_fish, form = ~Scientific_name), data = pruned_dat[pruned_dat$Class == "Fish", ], method = "ML")
summary(model18)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan + log_mass:log.lifespan 
>   Data: pruned_dat[pruned_dat$Class == "Fish", ] 
>        AIC      BIC    logLik
>   52.46258 58.75307 -21.23129
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                            Value Std.Error   t-value p-value
> (Intercept)            2.8425234 0.6822877  4.166165  0.0004
> log_mass              -0.0319071 0.1149743 -0.277515  0.7840
> log.lifespan           0.0608574 0.2242946  0.271328  0.7887
> log_mass:log.lifespan -0.0097255 0.0392377 -0.247860  0.8065
> 
>  Correlation: 
>                       (Intr) lg_mss lg.lfs
> log_mass              -0.420              
> log.lifespan          -0.568  0.445       
> log_mass:log.lifespan  0.421 -0.913 -0.694
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.52656470 -1.07403325 -0.87097331  0.01258059  0.90215128 
> 
> Residual standard error: 0.9651793 
> Degrees of freedom: 26 total; 22 residual
model19 <- update(model18, ~ .-log_mass:log.lifespan)
summary(model19)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass + log.lifespan 
>   Data: pruned_dat[pruned_dat$Class == "Fish", ] 
>        AIC      BIC    logLik
>   50.53509 55.56747 -21.26754
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                   Value Std.Error   t-value p-value
> (Intercept)   2.9137354 0.6060885  4.807442  0.0001
> log_mass     -0.0579381 0.0458255 -1.264320  0.2188
> log.lifespan  0.0222689 0.1581329  0.140824  0.8892
> 
>  Correlation: 
>              (Intr) lg_mss
> log_mass     -0.094       
> log.lifespan -0.423 -0.644
> 
> Standardized residuals:
>           Min            Q1           Med            Q3           Max 
> -1.5039949921 -1.0703688400 -0.8607380401 -0.0004501222  0.8154978887 
> 
> Residual standard error: 0.966526 
> Degrees of freedom: 26 total; 23 residual
model20 <- update(model19, ~ .-log.lifespan)
summary(model20)
> Generalized least squares fit by maximum likelihood
>   Model: new.log.TL ~ log_mass 
>   Data: pruned_dat[pruned_dat$Class == "Fish", ] 
>       AIC      BIC    logLik
>   48.5575 52.33179 -21.27875
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                 Value Std.Error   t-value p-value
> (Intercept)  2.949802 0.5379849  5.483057  0.0000
> log_mass    -0.053780 0.0343219 -1.566929  0.1302
> 
>  Correlation: 
>          (Intr)
> log_mass -0.529
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.50239115 -1.05995948 -0.87321945  0.01304271  0.83486397 
> 
> Residual standard error: 0.9669426 
> Degrees of freedom: 26 total; 24 residual
pruned_dat$Telomerase_activity[pruned_dat$Telomerase_activity == 0] <- "absent"
pruned_dat$Telomerase_activity[pruned_dat$Telomerase_activity == 1] <- "present"
pruned_dat$Telomerase_activity[is.na(pruned_dat$Telomerase_activity)] <- "N/A"

tel.act <- setNames(pruned_dat$Telomerase_activity, rownames(pruned_dat))
tel.act
>         Pygoscelis_adeliae         Tachymarptis_melba 
>                      "N/A"                      "N/A" 
>   Haliaeetus_leucocephalus       Phoenicopterus_ruber 
>                      "N/A"                      "N/A" 
>           Falco_sparverius      Ichthyaetus_audouinii 
>                      "N/A"                      "N/A" 
>            Hirundo_rustica           Branta_leucopsis 
>                      "N/A"                      "N/A" 
>             Cepphus_grylle           Rissa_tridactyla 
>                      "N/A"                      "N/A" 
>              Turdus_merula              Gallus_gallus 
>                      "N/A"                   "absent" 
>             Sterna_hirundo            Calidris_alpina 
>                  "present"                      "N/A" 
>      Haematopus_ostralegus    Aphelocoma_coerulescens 
>                      "N/A"                      "N/A" 
>        Chrysolophus_pictus              Fregata_minor 
>                      "N/A"                      "N/A" 
>                Parus_major                Gyps_fulvus 
>                      "N/A"                      "N/A" 
>        Strigops_habroptila      Oceanodroma_leucorhoa 
>                      "N/A"                  "present" 
>               Larus_fuscus      Tachycineta_albilinea 
>                      "N/A"                      "N/A" 
>        Colinus_virginianus         Fulmarus_glacialis 
>                      "N/A"                      "N/A" 
>          Macronectes_halli         Accipiter_gentilis 
>                      "N/A"                      "N/A" 
>          Amazona_amazonica          Buteo_jamaicensis 
>                      "N/A"                      "N/A" 
>            Vireo_olivaceus                  Sula_sula 
>                      "N/A"                      "N/A" 
>            Riparia_riparia  Passerculus_sandwichensis 
>                      "N/A"                      "N/A" 
>      Macronectes_giganteus            Buteo_swainsoni 
>                      "N/A"                      "N/A" 
>                Uria_lomvia     Aphelocoma_ultramarina 
>                      "N/A"                      "N/A" 
>        Tachycineta_bicolor           Diomedea_exulans 
>                  "present"                      "N/A" 
>           Lonchura_striata        Meleagris_gallopavo 
>                      "N/A"                      "N/A" 
>        Taeniopygia_guttata       Calonectris_diomedea 
>                  "present"                      "N/A" 
>  Acrocephalus_sechellensis          Anguilla_rostrata 
>                      "N/A"                  "present" 
>            Oryzias_latipes       Pseudocaranx_wrighti 
>                  "present"                      "N/A" 
>      Oreochromis_niloticus          Carassius_auratus 
>                      "N/A"                      "N/A" 
>            Cyprinus_carpio                Danio_rerio 
>                      "N/A"                  "present" 
>      Xiphophorus_maculatus       Xiphophorus_hellerii 
>                      "N/A"                      "N/A" 
>     Nothobranchius_furzeri    Nothobranchius_rachovii 
>                  "present"                      "N/A" 
>      Fundulus_heteroclitus               Gadus_morhua 
>                  "present"                  "present" 
>       Dicentrarchus_labrax   Acanthopagrus_schlegelii 
>                      "N/A"                      "N/A" 
>     Upeneichthys_vlamingii          Macquaria_ambigua 
>                      "N/A"                      "N/A" 
>  Lutjanus_argentimaculatus        Oncorhynchus_mykiss 
>                      "N/A"                  "present" 
>    Platycephalus_bassensis          Mustelus_asterias 
>                      "N/A"                      "N/A" 
>        Callorhinchus_milii   Urolophus_paucimaculatus 
>                      "N/A"                      "N/A" 
>          Squalus_acanthias          Torpedo_marmorata 
>                  "present"                      "N/A" 
>      Neoceratodus_forsteri            Setifer_setosus 
>                      "N/A"                  "present" 
>         Balaena_mysticetus      Eschrichtius_robustus 
>                   "absent"                   "absent" 
>     Giraffa_camelopardalis                 Ovis_aries 
>                   "absent"                   "absent" 
>          Rangifer_tarandus               Capra_hircus 
>                      "N/A"                      "N/A" 
>         Tursiops_truncatus        Camelus_dromedarius 
>                   "absent"                  "present" 
>                 Sus_scrofa          Muntiacus_reevesi 
>                  "present"                   "absent" 
>          Muntiacus_muntjak                 Bos_taurus 
>                   "absent"                   "absent" 
>     Zalophus_californianus            Crocuta_crocuta 
>                   "absent"                   "absent" 
>                Canis_lupus            Ursus_maritimus 
>                   "absent"                   "absent" 
>                Felis_catus            Ailurus_fulgens 
>                   "absent"                   "absent" 
>            Panthera_tigris                Meles_meles 
>                  "present"                      "N/A" 
>           Myotis_lucifugus      Tadarida_brasiliensis 
>                  "present"                  "present" 
>       Pteropus_rodricensis  Chaetophractus_vellerosus 
>                   "absent"                   "absent" 
>       Didelphis_virginiana       Atelerix_albiventris 
>                  "present"                  "present" 
>              Sorex_araneus          Procavia_capensis 
>                  "present"                   "absent" 
>         Lepus_californicus       Sylvilagus_aquaticus 
>                   "absent"                   "absent" 
>      Oryctolagus_cuniculus          Ochotona_princeps 
>                   "absent"                  "present" 
> Macroscelides_proboscideus       Galegeeska_rufescens 
>                  "present"                  "present" 
>        Ceratotherium_simum               Equus_grevyi 
>                   "absent"                   "absent" 
>             Equus_caballus            Tapirus_indicus 
>                   "absent"                   "absent" 
>    Myrmecophaga_tridactyla        Choloepus_hoffmanni 
>                   "absent"                   "absent" 
>               Homo_sapiens           Ateles_geoffroyi 
>                   "absent"                   "absent" 
>           Saimiri_sciureus               Pan_paniscus 
>                   "absent"                   "absent" 
>             Pongo_pygmaeus            Pan_troglodytes 
>                   "absent"                      "N/A" 
>            Gorilla_gorilla                Lemur_catta 
>                      "N/A"                   "absent" 
>          Macaca_nemestrina             Macaca_mulatta 
>                      "N/A"                   "absent" 
>        Macaca_fascicularis         Loxodonta_africana 
>                   "absent"                   "absent" 
>            Elephas_maximus          Castor_canadensis 
>                   "absent"                   "absent" 
>  Hydrochoerus_hydrochaeris       Sciurus_carolinensis 
>                   "absent"                  "present" 
>            Aplodontia_rufa      Heterocephalus_glaber 
>                   "absent"                  "present" 
>          Rattus_norvegicus               Mus_musculus 
>                  "present"                  "present" 
>                Tupaia_tana           Tupaia_belangeri 
>                  "present"                      "N/A" 
>         Alligator_sinensis Alligator_mississippiensis 
>                   "absent"                      "N/A" 
>           Zootoca_vivipara             Lacerta_agilis 
>                      "N/A"                      "N/A" 
>              Liasis_fuscus           Emys_orbicularis 
>                      "N/A"                      "N/A" 
>          Trachemys_scripta            Chrysemys_picta 
>                      "N/A"                  "present"
TA <- to.matrix(tel.act, unique(tel.act))
TA <- TA[pruned_data_tree$tip.label, ]

life.span <- setNames(pruned_dat$log.lifespan, rownames(pruned_dat))
log_mass <- setNames(pruned_dat$log_mass, rownames(pruned_dat))
telo.length <- setNames(pruned_dat$new.log.TL, rownames(pruned_dat))
plotTree(pruned_data_tree, ftype = "off", mar = c(3, 2, 2, 3))

tiplabels(pie = TA, piecol = c("white", "gray", "black"), cex = 0.22, offset = 4.3)

par(xpd = TRUE)

legend(x = 150, y = 0.2, legend = unique(tel.act), pch = 21, pt.bg = c("white", "gray", "black"), pt.cex = 1, bty = "n", title = "Telomerase activity", cex = 0.7, horiz = TRUE)


par(new = TRUE)
par(mar = c(3, 32, 2, 1.1))
barplot(life.span[pruned_data_tree$tip.label], horiz = TRUE, width = 1, space = 0,
        ylim = c(1, length(pruned_data_tree$tip.label))-0.5, names = "", las = 2, cex.axis = 0.5, axes = FALSE)

axis(1, at = round(seq(min(life.span), max(life.span), 1.5), 1), labels = FALSE)
text(round(seq(min(life.span), max(life.span), 1.5), 1), par("usr")[3] - 0.2, labels = round(seq(min(life.span), max(life.span), 1.5), 1), srt = 50, pos = 1, xpd = TRUE, cex = 0.5, offset = 1)
mtext("Lifespan \n (years)", side = 1, line = 1.6, cex = 0.5, font = 2)



Figure 3

## Ancestral state reconstruction of telomere size

fit <- fastAnc(pruned_data_tree, telo.length, vars = TRUE, CI = TRUE)

fit$CI[1, ]
> [1] 0.9673174 4.3349363
obj <- contMap(pruned_data_tree, telo.length, plot = FALSE)


##png("figure5.png", width = 7, height = 7, units = "in", res = 360)
##pdf("figure5.pdf")


plot(obj, ftype = "off", legend = FALSE, ylim = c(1-0.09*(Ntip(obj$tree)-1), Ntip(obj$tree)), mar = c(1, 0.1, 1, 5), lwd = 1.5)
add.color.bar(150, obj$cols,title = "Log telomere length (kb)", lims = obj$lims, digits = 3, prompt = FALSE, x = 0,
              y = 1-0.08*(Ntip(obj$tree)-1), lwd = 4, fsize = 0.6, subtitle = "")

tiplabels(pie = TA, piecol = c("white", "gray", "black"), cex = 0.16, offset = 5.7)
legend(x = 200, y = -5, legend = unique(tel.act), pch = 21, pt.bg = c("white", "gray", "black"), pt.cex = 1, bty = "n", title = "Telomerase activity", cex = 0.7, horiz = TRUE)

par(new = TRUE)
par(mar = c(3.6, 29.7, 3, 3))

barplot(life.span[pruned_data_tree$tip.label], horiz = TRUE, width = 1.07, space = 0,
        ylim = c(1, length(pruned_data_tree$tip.label))-0.5, names = "", las = 2, cex.axis = 0.5, axes = FALSE)

axis(1, at = round(seq(min(life.span), max(life.span), 1.5), 1), labels = FALSE)
text(round(seq(min(life.span), max(life.span), 1.5), 1), par("usr")[3] - 0.2, labels = round(seq(min(life.span), max(life.span), 1.5), 1), srt = 50, pos = 1, xpd = TRUE, cex = 0.5, offset = 1)
mtext("Log lifespan \n (years)", side = 1, line = 1.6, cex = 0.5, font = 2)

par(new = TRUE)
par(mar = c(3.6, 32.2, 3, 0.5))

barplot(log_mass[pruned_data_tree$tip.label], horiz = TRUE, width = 1.07, space = 0,
        ylim = c(1, length(pruned_data_tree$tip.label))-0.5, names = "", las = 2, cex.axis = 0.5, axes = FALSE)

axis(1, at = round(seq(min(log_mass), max(log_mass), 5), 1), labels = FALSE)
text(round(seq(min(log_mass), max(log_mass), 5), 1), par("usr")[3] - 0.2, labels = round(seq(min(log_mass), max(log_mass), 5), 1), srt = 50, pos = 1, xpd = TRUE, cex = 0.5, offset = 1)
mtext("Log mass \n (gr)", side = 1, line = 1.6, cex = 0.5, font = 2)

##dev.off()
## Correlated evolution under the threshold model


## Removing NAs from the dataset

pruned_dat$Telomerase_activity
>   [1] "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"    
>   [8] "N/A"     "N/A"     "N/A"     "N/A"     "absent"  "present" "N/A"    
>  [15] "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"    
>  [22] "present" "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"    
>  [29] "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"    
>  [36] "N/A"     "N/A"     "N/A"     "present" "N/A"     "N/A"     "N/A"    
>  [43] "present" "N/A"     "N/A"     "present" "present" "N/A"     "N/A"    
>  [50] "N/A"     "N/A"     "present" "N/A"     "N/A"     "present" "N/A"    
>  [57] "present" "present" "N/A"     "N/A"     "N/A"     "N/A"     "N/A"    
>  [64] "present" "N/A"     "N/A"     "N/A"     "N/A"     "present" "N/A"    
>  [71] "N/A"     "present" "absent"  "absent"  "absent"  "absent"  "N/A"    
>  [78] "N/A"     "absent"  "present" "present" "absent"  "absent"  "absent" 
>  [85] "absent"  "absent"  "absent"  "absent"  "absent"  "absent"  "present"
>  [92] "N/A"     "present" "present" "absent"  "absent"  "present" "present"
>  [99] "present" "absent"  "absent"  "absent"  "absent"  "present" "present"
> [106] "present" "absent"  "absent"  "absent"  "absent"  "absent"  "absent" 
> [113] "absent"  "absent"  "absent"  "absent"  "absent"  "N/A"     "N/A"    
> [120] "absent"  "N/A"     "absent"  "absent"  "absent"  "absent"  "absent" 
> [127] "absent"  "present" "absent"  "present" "present" "present" "present"
> [134] "N/A"     "absent"  "N/A"     "N/A"     "N/A"     "N/A"     "N/A"    
> [141] "N/A"     "present"
pruned_dat_not_NAs <- pruned_dat[!pruned_dat$Telomerase_activity == "N/A", ]
pruned_dat_not_NAs$Telomerase_activity
>  [1] "absent"  "present" "present" "present" "present" "present" "present"
>  [8] "present" "present" "present" "present" "present" "present" "present"
> [15] "absent"  "absent"  "absent"  "absent"  "absent"  "present" "present"
> [22] "absent"  "absent"  "absent"  "absent"  "absent"  "absent"  "absent" 
> [29] "absent"  "absent"  "present" "present" "present" "absent"  "absent" 
> [36] "present" "present" "present" "absent"  "absent"  "absent"  "absent" 
> [43] "present" "present" "present" "absent"  "absent"  "absent"  "absent" 
> [50] "absent"  "absent"  "absent"  "absent"  "absent"  "absent"  "absent" 
> [57] "absent"  "absent"  "absent"  "absent"  "absent"  "absent"  "absent" 
> [64] "present" "absent"  "present" "present" "present" "present" "absent" 
> [71] "present"
check2 <- name.check(pruned_data_tree, pruned_dat_not_NAs)
rm_phy2 <- check2$tree_not_data
rm_dat2 <- check2$data_not_tree
pruned_data_tree_not_NAs <- drop.tip(pruned_data_tree, rm_phy2)
pruned_dat_not_NAs <- subset(pruned_dat_not_NAs, subset = pruned_dat_not_NAs$Scientific_name %in% pruned_data_tree_not_NAs$tip, select = names(pruned_dat_not_NAs))

pruned_dat_not_NAs$Telomerase_activity <- as.factor(pruned_dat_not_NAs$Telomerase_activity)
str(pruned_dat_not_NAs)
> 'data.frame': 71 obs. of  20 variables:
>  $ Common.name               : chr  "Chicken " "common tern " "Leach storm petrel " "tree swallow " ...
>  $ Scientific_name           : chr  "Gallus_gallus" "Sterna_hirundo" "Oceanodroma_leucorhoa" "Tachycineta_bicolor" ...
>  $ Domain                    : chr  "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
>  $ Kingdom                   : chr  "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
>  $ Phylum                    : chr  "Chordata" "Chordata" "Chordata" "Chordata" ...
>  $ Class                     : chr  "Aves" "Aves" "Aves" "Aves" ...
>  $ Order                     : chr  "Galliformes" "Charadriiformes" "Procellariiformes" "Passeriformes" ...
>  $ Family                    : chr  "Phasianidae" "Sternidae" "Hydrobatidae" "Hirundinidae" ...
>  $ Genus                     : chr  "Gallus" "Sterna" "Oceanodroma" "Tachynicineta" ...
>  $ Species                   : chr  "gallus" "hirundo" "leucorhoa" "bicolor" ...
>  $ Endo_ectotherm            : chr  "endo" "endo" "endo" "endo" ...
>  $ Adult_mass_grams          : num  779.8 120 44.6 19 12 ...
>  $ log_mass                  : num  6.66 4.8 3.82 3 2.56 ...
>  $ Lifespan_years            : num  30 33 36 12.1 12 50 5 5.5 0.23 4 ...
>  $ Loglifespan               : num  1.48 1.52 1.56 1.08 1.08 ...
>  $ Average_Telomere_Length_kb: num  19.7 12.8 19 13.6 13.7 ...
>  $ logTL                     : num  1.29 1.11 1.28 1.13 1.14 ...
>  $ Telomerase_activity       : Factor w/ 2 levels "absent","present": 1 2 2 2 2 2 2 2 2 2 ...
>  $ log.lifespan              : num  3.43 3.53 3.61 2.57 2.56 ...
>  $ new.log.TL                : num  2.98 2.55 2.95 2.61 2.62 ...
head(pruned_dat_not_NAs)
>                               Common.name       Scientific_name    Domain
> Gallus_gallus                    Chicken          Gallus_gallus Eukaryota
> Sterna_hirundo               common tern         Sterna_hirundo Eukaryota
> Oceanodroma_leucorhoa Leach storm petrel  Oceanodroma_leucorhoa Eukaryota
> Tachycineta_bicolor         tree swallow    Tachycineta_bicolor Eukaryota
> Taeniopygia_guttata           Zebra finch   Taeniopygia_guttata Eukaryota
> Anguilla_rostrata            American eel     Anguilla_rostrata Eukaryota
>                       Kingdom   Phylum Class             Order       Family
> Gallus_gallus         Metazoa Chordata  Aves       Galliformes  Phasianidae
> Sterna_hirundo        Metazoa Chordata  Aves   Charadriiformes    Sternidae
> Oceanodroma_leucorhoa Metazoa Chordata  Aves Procellariiformes Hydrobatidae
> Tachycineta_bicolor   Metazoa Chordata  Aves     Passeriformes Hirundinidae
> Taeniopygia_guttata   Metazoa Chordata  Aves     Passeriformes  Estrildidae
> Anguilla_rostrata     Metazoa Chordata  Fish    Anguilliformes     Anguilla
>                               Genus   Species Endo_ectotherm Adult_mass_grams
> Gallus_gallus                Gallus    gallus           endo            779.8
> Sterna_hirundo               Sterna   hirundo           endo            120.0
> Oceanodroma_leucorhoa   Oceanodroma leucorhoa           endo             44.6
> Tachycineta_bicolor   Tachynicineta   bicolor           endo             19.0
> Taeniopygia_guttata     Taeniopygia   guttata           endo             12.0
> Anguilla_rostrata          Anguilla  rostrata           ecto           4032.0
>                       log_mass Lifespan_years Loglifespan
> Gallus_gallus         6.660319           30.0    1.477121
> Sterna_hirundo        4.795791           33.0    1.518514
> Oceanodroma_leucorhoa 3.819908           36.0    1.556303
> Tachycineta_bicolor   2.995732           12.1    1.082785
> Taeniopygia_guttata   2.564949           12.0    1.079181
> Anguilla_rostrata     8.302266           50.0    1.698970
>                       Average_Telomere_Length_kb    logTL Telomerase_activity
> Gallus_gallus                              19.70 1.294466              absent
> Sterna_hirundo                             12.80 1.107210             present
> Oceanodroma_leucorhoa                      19.03 1.279439             present
> Tachycineta_bicolor                        13.60 1.133539             present
> Taeniopygia_guttata                        13.70 1.136721             present
> Anguilla_rostrata                          12.50 1.096910             present
>                       log.lifespan new.log.TL
> Gallus_gallus             3.433987   2.980619
> Sterna_hirundo            3.526361   2.549445
> Oceanodroma_leucorhoa     3.610918   2.946017
> Tachycineta_bicolor       2.572612   2.610070
> Taeniopygia_guttata       2.564949   2.617396
> Anguilla_rostrata         3.931826   2.525729
pruned_data_tree_not_NAs
> 
> Phylogenetic tree with 71 tips and 70 internal nodes.
> 
> Tip labels:
>   Squalus_acanthias, Oncorhynchus_mykiss, Oryzias_latipes, Fundulus_heteroclitus, Nothobranchius_furzeri, Gadus_morhua, ...
> Node labels:
>   , '286', '79', '80', '87', '84', ...
> 
> Rooted; includes branch lengths.
name.check(pruned_data_tree_not_NAs, pruned_dat_not_NAs)
> [1] "OK"
names(pruned_dat_not_NAs)
>  [1] "Common.name"                "Scientific_name"           
>  [3] "Domain"                     "Kingdom"                   
>  [5] "Phylum"                     "Class"                     
>  [7] "Order"                      "Family"                    
>  [9] "Genus"                      "Species"                   
> [11] "Endo_ectotherm"             "Adult_mass_grams"          
> [13] "log_mass"                   "Lifespan_years"            
> [15] "Loglifespan"                "Average_Telomere_Length_kb"
> [17] "logTL"                      "Telomerase_activity"       
> [19] "log.lifespan"               "new.log.TL"
## Set the number of generations

ngen <- 5e6

## Run MCMC

mcmc.model <- threshBayes(pruned_data_tree_not_NAs, pruned_dat_not_NAs[ , c(17, 18)], type = c("cont", "disc"), ngen = ngen,
                          plot = FALSE, control = list(print.interval = 5e+05))
> Starting MCMC....
> generation: 500000; mean acceptance rate: 0.26
> generation: 1000000; mean acceptance rate: 0.24
> generation: 1500000; mean acceptance rate: 0.26
> generation: 2000000; mean acceptance rate: 0.22
> generation: 2500000; mean acceptance rate: 0.26
> generation: 3000000; mean acceptance rate: 0.24
> generation: 3500000; mean acceptance rate: 0.22
> generation: 4000000; mean acceptance rate: 0.24
> generation: 4500000; mean acceptance rate: 0.25
> generation: 5000000; mean acceptance rate: 0.22
> Done MCMC.
mcmc.model
> 
> Object of class "threshBayes" consisting of a matrix (L) of
> sampled liabilities for the tips of the tree & a second matrix
> (par) with the sample model parameters & correlation.
> 
> Mean correlation (r) from the posterior sample is: 0.478.
> 
> Ordination of discrete traits:
> 
>   Trait 2: absent <-> present
## Pull out the post burn-in sample and compute HPD

r.mcmc <- tail(mcmc.model$par$r, 0.8*nrow(mcmc.model$par))
class(r.mcmc) <- "mcmc"

hpd.r <- HPDinterval(r.mcmc)
hpd.r
>          lower     upper
> var1 0.1810489 0.7420291
> attr(,"Probability")
> [1] 0.9500012



Figure S1

## Profile plots from a Bayesian MCMC analysis of the threshold model

#png("figureS1.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figureS1.pdf")
plot(mcmc.model)

#dev.off()



Figure 6

## Plot posterior density

#png("figure6.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure6.pdf")


par(las = 1)
plot(density(mcmc.model), xlim = c(-1, 1.5))

## add whiskers to show HPD

h <- 0-par()$usr[3]
lines(x = hpd.r, y = rep(-h/2, 2))
lines(x = rep(hpd.r[1], 2), y = c(-0.3, -0.7)*h)
lines(x = rep(hpd.r[2], 2), y = c(-0.3, -0.7)*h)
box()

##dev.off()